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ABSTRACT 

We present theoretical analysis of the astrometric searches for extrasolar planets 
with the Space Interferometry Mission (SIM). Particularly, we derive a model for the 
future measurements with SIM and discuss the problem of reliable estimation of orbital 
elements of planets. For this purpose we propose a new method of data analysis and 
present a numerical test of its application on simulated SIM astrometric measurements of 
V Andromedae planetary system. We demonstrate that our approach allows successfull 
determination of its orbital elements. 

Subject headings: astrometry — methods: data analysis — planetary systems 



1. Introduction 

One of the most important and challenging goals of the Space Interferometry Mission (SIM, 
see http : / / sim . jpl . nasa . gov/) is astrometric detection of extrasolar planetary systems including 
Earth-like planets around stars from the solar neighborhood. High precision astrometry requires 
not only advanced technology but also adequately elaborated methods of data analysis. Among 
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others it is important to develop techniques allowing reliable detection of planetary signatures and 
extraction of the orbital elements. The aim of this paper is to discuss some of the problems related 
to this subject. Specifically, wc propose a method called the Frequency Decomposition (FD) to 
detect planets and help to obtain their orbital elements. This method has been successfully used 
for PSR 1257+12 timing observations (Konacki, Maciejewski & Wolszczan 1999) and 16 Cygni B 
radial velocity measurements (Konacki and Maciejewski 1999). Particular nature of the astrometric 
observations requires however some modifications of our original approach. In the paper we present 
a theoretical background of the method and an example of its application. 

The astrometric signal is a superposition of several effects of different magnitude and a proper 
analysis of the observations requires at least rough a priori knowledge of how different effects 
contribute to the signal. However, these effects depend on the parameters (such as number of 
planets, their eccentricities) which are unknown in advance. So what we propose is a two-step 
analysis (1) FD to understand the basic properties of the signal (i.e determine the number of 
planets and approximate values of their orbital elements) and (2) least-squares fit based on a 
proper model and the starting values of parameters derived from the previous step to refine the 
parameters and obtain their uncertainties. 

The basic idea of FD is the following. With few exception (proper motion, long period planets) 
the processes contributing to the signal are periodic. Therefore our astrometric signal can be 
successfully modeled as a multiple Fourier series plus a polynomial of certain degree (to account for 
the proper motion and long period planets). FD is a numerical algorithm to obtain the estimates 
of frequencies, amplitudes and phases of such model (Konacki, Maciejewski k, Wolszczan 1999). 
Let us note that such approach, contrary to the usual least-squares method, allows us to analyze 
the data without assuming any physical model. Subsequently we interpret derived parameters. 
We decide how many planets are present in the system and calculate their orbital parameters (as 
one can derive analytical formula expressing amplitudes and phases as functions of the orbital 
elements). This is especially useful for multiple planetary systems where deciphering the number of 
planets may be tricky (e.g. two planets in circular 2:1 resonant orbits may mimic one planet in an 
eccentric orbit, see Konacki and Maciejewski 1999). Our approach can be also helpful while trying 
to determine whether we observe an astrometric displacement from a planet in 1-yr orbit or anmial 
parallax since the parallactic motion has its own specific Fourier expansion constrained by SIM 
orbit. Finally, we can use these findings to perform the 'traditional' least-squares fit. We believe 
that such approach allows us to make more justified hypothesis about the data and in consequence 
lead to reliable results. 

The plan of our paper is the following. In section 2 we derive a detailed model of the SIM 
measurements. In section 3 we investigate Fourier properties of the orbital motion. In section 4 we 
discuss our approach to the analysis of SIM data and finally in section 5 we perform some numerical 
tests to show how our method works in practice. 
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2. Modeling delays 

SIM measures relative positions of stars using Michelson interferometers. A single measurement 
with SIM gives the projection of direction to the star s onto the interferometer baseline vector B. 
The measured quantity is the optical pathlength delay between the two arms of the interferometer 
(Shao & Baron 1999) 

d = B -s + c + e, (1) 

where c is the zero point of the metrology gauge and e represents measurement uncertainty. The 
search for extrasolar planet is performed in so-called narrow angle mode where delays toward 
two stars (called target and reference) within 1° are measured and compared. For this kind of 
observation the measured quantity becomes the relative delay 

= B • (si - S2) + e, (2) 

where si and S2 are directions to the target and reference stars, respectively. Such narrow-angle 
measurement gives the angular separation between the stars and offers higher accuracy as many 
errors scale with the angular distance. 

The direction to a star S = S{t) from the Solar System Barycenter (SSB) is changing with 
time due to the proper motion and presence of companions. These two effects, we model in the 
following way 

S(t) = So + <5S^(t) + <5Se(t), (3) 

where So is the direction toward the star at epoch to; '^S^j(t) and 5Sc{t) describe changes due to the 
proper and orbital motion, respectively. In order to properly calculate these changes let us assume 
that the SSB radius vector of each star is given by 

R = Ro + 5R, where ||(5R|| < ||Ro||, (4) 

then up to the first order in ||(5R||/||Ro || 

S = ^ = So + (5SW, (5) 



where 



||Roll llR-i 
The correction 5S^^^ can be written in the form 



^" 5S(^) = -^^Sq X (So X SK) , (6) 



<5S(i) = -— [<5R - So(So • <5R)] = j^5R±, (7) 
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which means that it only depends on the component of STL perpendicular to Sq. It turns out that 
sometimes the first order correction is not sufHcient. Therefore we need to derive and analyze also 
the second order term SS^'^^ given by 



3 / SK Y 1 /||(5R"^ ^' 
t: ^0 



2 V / 2 V Ro 



"^'S..^! (8) 



Roll \ IIR-ol 



If we represent (5R as a sum of two components, perpendicular and parallel to Sq, 

6R = 6R± + 6R\i (9) 

can be written as 



(10) 



As we show in section 2.3, dS^'^^ is especially significant for nearby stars with large proper motions. 
For such stars we have 

5K±=\'Tt (5R||=Vi?i (11) 

where V^- and Vij arc respectively transverse and radial velocity of the star. Thus if the star has 
a significant proper motion, through astrometric observations we can detect angular displacement 
56 (second term in equation (10)) 

^"-r",," "ifriNsi^'^ (12) 

||Ro|| IjRoll Ijl^oll Ijl^oll 

due to the radial velocity. This effect is called perspective acceleration. The other term from the 
equation (10) has an interesting property. Namely, it can be shown that it does not change the 
angle between S{t) = Sq + 5S^^)(i) + SS^'^\t) and Sq (i.e. current and initial position of the star). 
In other words, if we had a direct way to measure this angle, we would not observe any effect from 
that term. However, since we measure all angles through the equation (1) and model the unit vector 
toward the star with S{t) = Sq + (5S(i)(i) + 5S^'^\t), the term -i (||(5R_l||/||Ro||)^ Sq is necessary. 
Specifically it affects the length of the vector S and helps to keep it normalized within the accuracy 
of the second order approximation. Further details concerning second order corrections we discuss 
in section 2.3. 



2.1. Local frame and baseline vector orientations 

In order to obtain explicit form of the delay d we need to calculate a scalar product B • S. 
The value of this product does not depend on a chosen reference frame. Thus, depending on our 
needs we can express vectors on the right hand side of (3) in different ways. It is convenient to 
introduce a local right hand orthonormal frame at the point Sq on the celestial sphere (see Fig. 1). 
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This frame is connected with the classical equatorial spherical coordinates and is defined by the 
unit vectors {60,65,6^}. In SSB equatorial frame coordinates of these vectors are the following 

= (— sinQ:,cosa, 0), = (—sin (5 cos a, — sin (5 sin a, cos (5), (13) 

Br = So = (cos 5 cos a, cos (5 sin a, sin (5), (14) 

where {a, 5) = {ao, 5o) is the right ascension and declination of a star at to. 

One can determine the relative position of the target and reference star using two interfer- 
ometers or, as it is planned for SIM, by performing two measurements with one interferometer for 
two non parallel orientations of its baseline, Bj, i = 1,2. For each orientation, the baseline can 
be represented as a sum of two vectors, B^',B^, parallel and perpendicular to the initial direction 
toward the target star, Sq. Since we have S{t) = So + dS{t) where SS{t) is a displacement tangent 
to So, the delay can be written as 

d = Bi- S{t) + c + e = B;' • So + B,^ • 5S{t) + c + e = do + Ad(i) + c + e (15) 

where do = B'.' • Sq and Ad{t) = B^ • 6S{t). Clearly, from the planet detection point of view the 
important term is Ad{t). Assuming that the measurement uncertainty, e, is independent on the 
baseline orientation the most favorable orientation of Bj is Bj = B^. For such orientation Ad{t) 
is the largest possible for a given length of the baseline. Moreover, the baseline orientations should 
be perpendicular. This way the the covariance ellipse on the sky will be circular and there will not 
be any direction on the plane tangent at So in which the measurements are more accurate than 
in others. Therefore, for all further considerations we assume that all observations are made with 
two orthogonal and fixed baseline orientations Bi and B2 which are perpendicular to the initial 
direction toward the target star. Additionally, to simplify the equations we assume that Bi is 
parallel to Ba and B2 is parallel to eg. 



2.2. Proper motion, parallax and companions 

The proper motion is a projection on the sky of the motion of a star with the velocity V 
and within the first order approximation astrometrically only its transverse component Vt = 
VaGa + Vses is observable. Thus using simple arguments we find that 

6Sn{t) = ■K{Vatea + Vstes) = cos 5 /j^atea + M<5*e5, (16) 

where 

Va = Y-ea, Vs = \'-es and i^a = —{to), fis = —{to). (17) 
and TT = 1/D* where D* = ||Ro|| is the SSB distance to the star. 
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If the star has companions then the proper motion refers to the motion of the mass center of 
the system and the first order correction due to the orbital is given by the fohowing equation 

SScit) = TT [K(i)ea + RKt)es] , (18) 

where R* = (i?* , R^, R*) denotes the radius vector of the star with respect to the barycenter of its 
system in the local frame {ea,es,er}. 

If the interferometer is located at Ro(0 in SSB frame then the observed direction toward the 
star is 

s(t) = S(t) + n(t), (19) 

where Il(t) is the parallactic displacement 

n(t) = 7rSo X (So X Ro(t)) (20) 

obtained from the equation (6) by substituting — Rq for (5R. In the local frame the parallactic 
displacement can be written in the following form 

n{t) = n[na{t)ea + Ils{t)es], (21) 

where 

n«(t) = -Ro(t) • ea = Xo{t)sma - Yo{t)cosa, (22) 

^6 (t) = — Ro {t) ■ eg = Xq (t) sin 5 cos a + Iq (t) sin S sin a — Zq (t) cos S (23) 

and {Xo{t),Yo{t), Zo{t)) are the coordinates of SSB vector Ro(i)- The expressions for 8^(1), Sc{t) 
and n(t) come directly from the formulae (6), (7) and thus represent a first order approximation 
with respect to vr. 

Now, using (3)-(20) and assuming that the measurements have been already corrected for 
aberration and gravitational lensing we can rewrite the delay equation (1) in the form 

d = (f + d'^t + d'^ ■ Ro{t) + ■ R*(t) + c + e, (24) 

where 

d° = B ■ So, d^" = ttB ■ Vt = h ■ iJ,, fi = {/J-a cos 6, us), (25) 

h = {B-ea,B-es), d'^ = -tt [B - (B • So)So] , (26) 

d^ = 7rb, R*(t) = (i?*(0,i?f(t)). (27) 

Using (24) we obtain similar formula for the relative delay 

D = D° + D^'t + 'D'' ■ Ko{t) + d^ • R*(t) + e, (28) 
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where 



D'^ = [7ri(B ■ Sj)Sj - 7r2(B ■ S§)S§] - (tti - 7r2)B. 



(29) 



where the indices 1,2 refer to the target and reference star respectively. In the above we assumed 
that only the target star has companions so as d'^ refers to the local reference frame of the target 
star. Formula (28) plays the fundamental role in our consideration. It explicitly shows the structure 
of the observed signal that consists of a dominant linear trend modulated by periodicities due to 
the motion of the interferometer and the companions of the target star. 

According to our assumptions a single observation is done for two orthogonal baseline orien- 
tations Bi and B2. Thus a single observations is given as a two component vector D = {Di,D2) 
of the relative delays. Each component Di has the form (28) where coefficients D^, D!^, DJ and d? 
are calculated with the formulae (29) and B = Bj, i = 1,2 respectively. 



The above considerations represent first order approximation which is sufficient for most as- 
trometric measurements. However SIM is expected to deliver unprecedented 1/zas precision in the 
narrow-angle mode and thus it is important to understand limitation of the model (29). It can be 
accomplished by analyzing higher order terms. For the baselines perpendicular to So, the second 
order corrections correspond to the term that is responsible for the actual angular displacement 
(see equation (12)) and we have 



They can be calculated if we put SR = -Ro{t) + R*(i) + Ry (i) where Ry (t) = Vr t + Vi? t. We 



2.3. 



Second order corrections 




(30) 



obtain 



-^<5S(2) = -||R^(t)||R^(t) + ||R:(i)||Rl(t) + ||R^(i)||R^(t)+ 

+ (llR^(i)IIRl(t) - ||R:(t)||Ri(i)) + (llR^(t)llR^(i) - ||R^(t)||Ri(t)) + (31) 
+ (||R:(t)||R^(t) + ||R^(t)||Ri(t)) 



and 



||5S(2)|| < 7r2(||Ro(t)f + \\K*{t)f + \\Kv{t)f+ 



(32) 



2Ro(t) •R''(t) -2Ro(t) •Ry(i) + 2R*(t) •Ry(t)) 
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As we can see there are two types of second order corrections. The first type includes the second 
order corrections due to the proper motion, parallax and companions. For the proper motion it is 
easy to calculate its exact value 



= T = ^VtVr AT^ (33) 

4 llRoil ||Ro|l 4 ^ « ^ ' 



where AT is the time span of the mission (and we assumed that to is at half of AT), ||Vx|| = Vt, 
W^rW = III order to learn about the magnitude of this correction, we calculated its value for 
the sample of 150 stars from the Hipparcos catalogue with the largest proper motion (see the Inter- 
net location http : / / astro . estec . esa.nl/SA-general/Pro jects/Hipparcos/hipparcos .html). 
The results are shown in Fig. 3. As one can see AS^ is indeed significant for such stars and without 
any doubts has to be included into the model. For the remaining corrections due to the motion of 
the interferometer (i.e. the second order parallactic correction) and the presence of a companion 
we have the following upper limits 

7r2||R^(t)||||R^(t)|| < 7r2||Ro(t)f < AS^ = ttVo, 

(34) 

7r2||R;;(t)||i|Rl(t)|| < 7:^\\K*{t)f < AS, = tt^ (a./(l - e,)f 

where ao is SIM's orbit semi-major axis, a* semi-major axis of the star's orbit and its eccentricity. 
If we assume that SIM semi-major axis is 1 AU then 

ASc ~ 5 X 10~ -P;^ aas for planetary companions 



ml, P^i^ (35) 
ASc ~ 4.9 -p: jias for stellar companions 

dicM^,^ + niMjMM^r/Hl - e)2 

AStj- « 4.9 jias 

where dpc is the distance to the star in parsecs, Mmq mass of the star in solar masses, mjvfg mass 
of the companion in solar masses, iriMjup mass of the companion in Jupiter masses, Pyr orbital 
period in years and e eccentricity. 

The other type of second order corrections includes all " mixed" terms 

7r2||||R;;)(i)||Rl(t) - ||R:(t)||R^(t)|| < 27r2 \Ro(t) ■ R*(i)| < AU, = 27r''aoaJ{l - e^), 

TT^WKimKit) + IIRy WIIRiWII < 27r2|R*(t) • Ry(t)| < A*e = 7rVAra*/(l - e^, (36) 

7r2||||R^(t)||R^(t) - ||R^(t)||R^(t)|| < 27r' \Ko{t) ■ Kv{t)\ < AU^ = ir^aoVAT 

where V = (V^ + ^)^''^- The value of AII^ has been calculated for the same sample of stars as 
in Fig. 3. Again, one can observe that AII^ is significant and a proper model has to take it into 
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account (see Fig. 4). We also find that 

m P^^^ 

AIIc 9 X 10~^ '!}'!^ — lias for planetary companions 

dl,M%l{l-e) 

rriM. P^^^ 

AIIc ~ 9.7 r-rr £qj, stellar companions 

■muQlMuofl'^i}- - e) 



(37) 



and 



A^'c !=a 0.1 ViooATyr /las for planetary companions 



p2/3 

A^'c ~ 102.3 ViooATyr —r- @__yji ^ £qj, stellar companions 

d2,M2^(l + mMe/MMj2/3(i_e) 



(38) 



where Fioo is the velocity of star in hundreds of km/s and ATyr is the time span of the mission in 
years. 

Finally, let us shortly discuss the magnitude of third order terms. Obviously, they will be 
detectable only for SH due to the proper motion. Thus we can estimate that the resulting angular 
displacement 6S^^^ is 

\||Ro||/ 8 8 

Its value for the sample of stars from Figs. 3-4 is presented in Fig. 5. As one can see in few cases 
this third order term can be larger than Ifias. 

The above analysis clearly shows that a variety of second order effects and possibly in few 
cases third order effects will be detectable with SIM. Although throughout the rest of this paper 
we use only the first order model presented in sections 2.1-2 to simplify our considerations, in real 
applications a correct model must include higher order effects. They can be easily derived given 
the theoretical background presented in section 2. 



3. Orbital motion 

Let us assume that the motion of N planets and their star is described in the barycentric 
system. From the definition of such system, we have the following relation for the radius vector of 
the parent star 

1 ^ 

^* = -iw-^"'^'^^' ^^^^ 
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where Hj are radius vectors of planets and M*,, rrij are the mass of the star and the j-th planet, 
respectively. In the first approximation, the motion of planets can be described by means of the 
following equations 

^-^/^^-^E^' J = h---,N, (41) 



where 



"J/ 

and the motion of the star can be obtained from the equation (40) . 



3.1. Elliptic motion and its expansion 

Solutions Rj(i) of the equations (41) belong to the family of Keplerian orbits among which 
elliptic orbits are of particular interest for farther analysis. Therefore, let us remind their basic 
properties. 



The radius vector Hj = R(t) of a planet moving in an elliptic orbit is given by 



where 



R(i) = P a(cos E(t) - e) + Q a^l - sin E{t), 



P = 1 cosa; + m sina;, Q = — 1 sina; + mcosw. 



(42) 





cos Q 




— cos i sin $7 


1 = 


sin 0, 


m = 


cos i cos Q, 









sini 



The eccentric anomaly E = E(t) is an implicit function of time through the Kepler equation 

E -esinE = M, (43) 

where M is the mean anomaly 

M = n{t-Tp), n=—, (44) 



and P is the orbital period of a planet. The remaining parameters a,e,ijj,^l,Tp are the standard 
Keplerian elements — semi-major axis, eccentricity, longitude of pericenter, longitude of ascending 
node and time of pericenter. 
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The functions cosE' and sin£^ are periodic with respect to M. and can be expanded in the 
Fourier series 

cos£' = — ^ ^Jk-i{ke) cos{kAi) , 



sin£' = ^ ^Jk^i{ke)sm{kM), 



k 

keZa 

where Jn{z) is a Bessel function of the first kind of order n and argument z; Zq denotes the set of 
all positive and negative integers excluding zero, and e G [0, 1). Thus, using the equations (42) and 
(45), we obtain 

R(i) = R° + A ^ ^Jk-i{ke) cos(fcM) + B ^ ^Jk-i{ke) sin(fcM), (46) 
keZo keZo 

where 

1 3 
R(t) = -R, R° = — P e, A = P B = Q Vl - e^. 

a 2 

It can be written in the following complex form 

%) = R°+ ^ ©fee''^^, (47) 
keZo 

where 

&k = ^{F-ik,e)A-iF+{k,e)B), (48) 

and 

F± {k,e) = Jk-i (ke) ± Jk+i (ke) . (49) 

Eventually, using (44) and (47) , we obtain the Fourier expansion of R(^) 

R(t) = R° + ^ Afe e'*^"*, where A^ = ©^6-^*^^*^^ (50) 
keZo 

Let us define the following quantity 

Ai = ^l^fl^L for / = 1,2,3, and A; > 0, (51) 

i.e. the ratio of amplitudes of two successive harmonics, where A* is the i-th component of vector 
Aj . From the properties of Bessel functions we have 



Aiie) ^ 



k + 1' 



\ e^Aiy[4{ke)f + {Biy[Jk{ke)Y 



(52) 



where J^(^) indicates the derivative of a Bessel function Jn{z) with respect to z. It can be proved 
that for all e G (0, 1), I G {1,2,3} and A; > we have .4|j(e) < 1. It means that the expansion 
of R(i) has an important property — moduli of successive harmonics of each of coordinates of R(t) 
decrease strictly monotonically with k. 
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3.2. Real expansion 



Given the equations from the previous section, it is possible to derive the real expansion for 
every component of the vector R(i) = {Ri{t), R2{t), R3{t)). Namely, we can express this vector in 
the form 

oo 

R(t) = rO + ^ (c'' cos{knt) + S*^ sin(feni)) (53) 
fe=i 

which is more convenient in numerical applications. Using (47), (48) and (50) we find 

1 



c = 



k I 
1 

k 



PF^{k, e) cos{knTp) - Q^/l - e'^F+{k, e) sin(A;nTp) 



PF_(fc, e) sin(A;nTp) + Qy/l - e^F+{k, e) cos{knTp) 



(54) 



From the above formulae immediately follows that amplitudes of successive harmonics are given by 



{Dff = {C^f + {S^f = 1 [P^F^{k,ef + QUi - e')F+{k,ef] , Z = 1,2,3. 



(55) 



In applications it is convenient to have these expressions in an explicit form 
1 

P L 
1 

1 



(^Di f = ^ F_(fe,e)^ (l -sin^isin^ri) + F(A;,e) (cosfisina; + coszsinJ^cosw)^ 
(I?2)^ = F-{k, e)^ (l — sm^ icos^ + F(A;,e) (sinrisina; — coszcosficoso;)^ 



(56) 



where 



, „ \F_ (k, ef + F{k, e) cos^ uA sin^ i, 

F{k,e) = {l-e^)F+{k,ef -F^{k,ef. 



3.3. Approximate formulae for small and moderate eccentricities 

Since small and moderate eccentricities are more probable it is useful to have approximations 
of the expressions from the previous section. Namely, using known expansions for Bessel functions 
we obtain the following formulae 

^*(*'') = Jkhy. (y)"" +0(e'«). ^'F,(k,e) = ^ (I)"" +0(e'«) (57) 
Subsequently 



k{k-\)\ V 2 



ke 



[Icoswfc + msinu;jk] + 0{e^^^) 



1 fke\'-' ^''^ 
}^(^j^_-^y_ [-Isi^'^fe + mcoswfc] + 0{e''+^). 
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where cok = co — knTp. Finally, we obtain the expansions for the amplitudes 



ke 



k{k- 1)1 V 2 



1 



ke 



k{k-l)\ \ 2 



1 



ke 



k{k-iy. \ 2 



k-l 



k-l 



(1 -sin^zsin^ n)+0{e'^''), 
(1 -sin^zcos^J^) +0(6^'=), 



(59) 



sm^i + Oie'"). 



2k\ 



Prom the above we can obtain the harmonic expansion for a circular orbit. Namely we find that 
R° = and C^^ = S'^ = for A; > 1. While for A; = 1 



C'^ = [1 cos a; + m sin a)] , 
S'^ = [— Isina; + mcosa)] . 



(60) 



where ui = —nTp. 



These equations will be especially useful for deriving orbital elements from the coefficients 
C'^jS'' obtained through the analysis of observations. We discuss this issue in the next section. 



4. Data analysis 

For the tests we assume the following SIM observing scenario. At the moments tj and tj+i the 
relative delay between the target and reference star is measured for two orthogonal baseline orienta- 
tions Bi and B2. Such measurement gives a two dimensional delay vector Dj = (Di(ij), D2(ii+i)) 
and is repeated N times over the time span of the mission, AT. As a result we obtain a two di- 
mensional time series V = {Di,i = 1, . . . , A^}. The goal of the data analysis is to detect planetary 
signatures in V and derive the orbital parameters of planets. We solve this problem in two steps. 
First wc perform Frequency Decomposition (FD) of the time scries T>. The aim of this step is 
to understand the basic properties of V i.e. determine the number of planets and estimate their 
orbital parameters. The second step is the least-squares analysis based on a specific physical model 
established in the previous step. Its is aim is to obtain accurate values of the orbital elements and 
their uncertainties. These two steps are described in the following sections. 
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4.1. Harmonic model 



From the theoretical considerations of sections 2 and 3 it foUows that relative delays can be 
modeled by means of the following expression 



D = D° + D'^ i + ^ [C'^'*^ cos{nokt) + S'^'^ cos{nokt) 

k=l 

N oo 

+ ['^^''^ cos(njfei) + S^'^ cos{njkt) 



+ 



(61) 



j=i k=i 



where A'' denotes the number of planets, no and nj denote the mean motion of SIM and j-th planet, 
respectively. Such equation comes directly from the fact that the motion of the interferometer and 
the motion of planets can be expanded into the Fourier series. Consequently, the above formula is 
used to describe V and special numerical algorithm is used to obtain the parameters 
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(62) 
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no, 



j = l,...,N, k = l,...,Kj, 



This algorithm has been described in great detail in Konacki, Maciejewski &: Wolszczan (1999). 
Let us only note here that in practice, due to limited accuracy of measurements and the fact that 
the amplitudes of subsequent harmonics decrease monotonically, the expansions of (61) are finite 
and the (finite) number of harmonics Kj depends mainly on orbital eccentricities and measurement 
errors. Using our algorithm wc can determine the number of planets N, the number of detectable 
harmonics Kj and determine the basic frequencies and coefficients of (61). 

In fact, we can assume that the mean motion of SIM no as well as the other elements of its 
orbit are known. In other words Ro(i) is known and we can use the following more constrained 
version of the formula (61) 



D = D" + D^t + 



N Kj 

Ro(0 + XI [^^''' cos{njkt) + S^'^ cos{njkt) 
j=i k=i 

This way instead of several parameters C^'^, S'^'^, no we have six parameters since 



r)7r r)7r r)7r 
-^115-^125 -^13 
7-)7r 7-)7r 7-)7r 
-^^215 -^^225 ^23 



(63) 



(64) 



In order to have a better understanding of the parameters of (63) let us express them by means 
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of the quantities introduced in section 2. Using (28), (29) and (50) we find that 
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^2 



where D^, and DJ^ are the quantities defined by (29) and calculated for B = Bj, i = 1,2; H^'^ 
denotes R*^ in the expansion (53) for the orbit of the j-th planet and aj is the semi-major axis of 
the j-th planet. The coordinates of R-' = {R\, -R3) are expressed in the local frame {ca, e^, e^}. 
This way 



d'^ • R^' = 7ri(Bi • ea)R{ = 7riBiR{, for Bi 
d'=-R^ =7ri(B2-e5)^^ = 7riS2^^, for B2 



(66) 



since Bi = 5i e^, B2 = B2 es where Bi is the length of the baseline vector Bj. Similarly we have 



,k 



BiCi 
B^Ci''' 



mi 



,k 



BiSi 
B^Si''' 



(67) 



where C^^ = (C^, C7^■'^ C^''') and = {Si''' , S^''' , S^''') are C^^ and S*^ coefHcients of expansion 
(53) for j-th planet expressed in the local frame. 



Now let us assume that after performing FD, we obtained the parameters 

D°, D^^, D^, c^■•^ s^■•^ j = i,... ,n, k = i,... ,Kj 



(68) 



where A'^ is the number of planets (i.e. the number of basic frequencies detected) and Kj is the 
number of detected harmonics for each planet. The first question is if we can derive the canonical 
parameters like a, 5, tt for the target and reference star from D'', D^, D'^. Unfortunately this 

is not possible, at least without additional assumptions. Obviously it is a direct consequence of 
the relative measurements we perform. Thus we can only derive (Sj — Sq) as well as differential 
proper motion and differential parallax. In fact it is possible to chose such a reference star that the 
differential parallactic displacement has an amplitude close to zero. It suffice to have a reference 
star with the parallax similar to the parallax of the target star since by assumption these two stars 
are close to each other and their parallactic displacement is very similar. This way we can remove 
a strong parallactic component from our observations. On the other hand we do not need the exact 
values of the canonical parameters. We only have to properly remove the respective effects in order 
to be able to detect putative planets. 

The remaining question is if we can derive the orbital elements from C''"'^, S^'*^. This task is 
relatively simple. Namely, given that we have detected at least two terms (basic frequency and its 
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first harmonics), all orbital elements can be derived from the equations of section 3.3. For planets 
with only the basic frequency detectable, we assume a circular orbit and then the other elements 
can also be found. This procedure we demonstrate in section 5. 

4.2. Standard model 

The harmonic model allows us to describe the data without a priori knowledge of the target star 
parameters and its planetary system. In the same time it allows to derive all important information 
— especially the number of planets and estimates of their orbital elements. With such knowledge 
we are ready to perform the standard least-squares analysis in which we must specify the model 
and supply good initial conditions for the fit. 

The standard model has the following form 

N 

T> = & + f>i't + W ■ Ro(t) + 

where Ri , R2 are coordinates of the Keplerian motion vector R given by the equation (46) . The 
parameters of such model are 

DO, D^, %,rp,,-,ej,i,-,cc;,-,J],-,P,-, j = l,...,N, (70) 

where Tpj,ej,ij,LOj,^}j are the Keplerian elements of the j-th. planet, Pj is its orbital period and 
the parameter aj is defined in the following way 

^ rrij 
5. Numerical tests 

For the tests we chose v And with its two outer planets (Butler et al. 1999). All real and 
assumed astrometric and orbital parameters are presented in Table 1. We also found a reference 
star HD 10032 which is about 0.°7 away from v And. Its astrometric parameters are in Table 2. SIM 
is assumed to move in an orbit similar to the orbit of the Earth (see Table 3). For these stars we 
simulated N = 200 measurements of relative delays {Di,D2) (for two baseline vector orientations 
Bi and B2) randomly distributed over the time span of 10 years. In both cases the length of the 
baseline vector was 10 meters and a measurement error with cr « 50 pm was assumed (i.e. l/ias 
in angular displacement). Since by assumption Bi is parallel to e^ and B2 to e^, the delay Di 
corresponds to an angular distance between v And and HD 10032 in right ascension and D2 to an 
angular distance in declination. 



aj Bi Ri (t, T^,j , Cj , ij , uj , , Pj] 
Gj B2 R2{t, Tpj,ej,ij,(jjj, Pj] 



(69) 



(71) 
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5.1. Second order effects 

Before we proceed with the analysis it is interesting to discuss second order effects present 
in simulated observations. Since v And is a nearby star with large proper motion we can expect 
significant contribution from this star (our reference star HD 10032 is quite distant and thus all 
second order effects are mainly due to v And). One can analyze these effects by means of the 
formulae from section 2.3 or simply apply the standard model (69) with the parameters precisely 
computed from assumed parameters of v And, HD 10032 and SIM (Tables 1-3) and examine the 
resulting residuals. This procedure gives the residuals presented in Fig. 6. As we can see the 
second order effects are dominated by a variation quadratic in time (Fig. 6 a,b) . This effect is due 
to perspective acceleration 

^IIV^IIVtAT^ (72) 

thus if we assume that the radial velocity, Vr, of v And and HD 10032 is zero it will disappear 
and reveal another second order effect of smaller magnitude (see Fig. 6 c,d). This effect is due to 
the following term 

7r2||R^(t)||R^(i)-7r2||R^(t)||R^(i) (73) 

i.e. the mixed term of parallax and motion of the star. 

Finally let us note that if we allow the parameters of the model (69) to vary, as usual during 
the process of least-squares fit, this first order model will try to minimize the residuals as presented 
in Fig. 6 e,f. 

5.2. PVequency Decomposition and standard model 

First step in our analysis of sinuilatcd relative delay measurements is the Frequency Decom- 
position (FD). Here we model the data with the less constrained formula (61) to show how the 
parallactic motion contributes to the data. From assumed parameters of the stars and SIM we can 
compute amplitudes of basic terms and their harmonics. They are shown in Fig. 7. The main idea 
of FD is subsequent removal of effects with decreasing magnitudes (for all details of the method 
see Konacki, Maciejewski & Wolszczan (1999)). This process is demonstrated in Fig. 8 and 9 for 
Di (i.e. for delays measured with the baseline vector orientation Bi). As one can see the most 
significant part of delay variations comes from the proper motion of both stars (Fig. 8a), then we 
can detect the basic term of the parallactic motion (Fig. 8b), the basic term of the planet II (Fig. 
8c), first harmonic of the parallactic motion (Fig. 8d), the basic term of the planet I (Fig. 9e), 
first harmonic of the planet II (Fig. 9f), second harmonic of the planet II (Fig. 9g) and finally first 
harmonic of the planet I (Fig. 9h). The values of respective parameters S-''^, C^'^ are presented in 
Table 4. They are sufficient to derive initial estimates of the orbital elements of planets I and II. 
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Namely, from the approximate equations of section 3.3 we can find that 



^ ^2 /(5l!)!±i3:!)! 



„„„ on „;„2 , 



^1 9 9 • — ~ - 9 



(74) 



— — \ — = a, sm ^Ij cos fi,- sm^ ij 

and, together with analogous formulae for first harmonics, easily determine the orbital elements. 
They are show in Table 5. As one can see this procedure gives quite accurate values of the orbital 
elements. However we use them only as initial values for the least-squares fit with the standard 
model (69) to obtain the final parameters presented in Table 6. 



5.3. Conclusions 



The above test demonstrates that our approach allows us to determine the orbital elements 

with high confidence, at least in this particular case. It is interesting to note that with FD we arc 
able to estimate the orbital elements without using the entire information present in the simulated 
data set (the residuals from Fig. 9h are well above the assumed measurement error). Surprisingly 
this estimation is quite accurate and as demonstrated is perfectly sufficient as an initial guess of 
the parameters for the standard least-squares analysis. This is a very promising result since the 
difficult problem of good initial condition is usually solved by means of quasi-global techniques 
which are very demanding numerically and still may lead to unreliable results. Thus we believe 
that our approach constitutes safe and efficient solution to the problem of planets detection with 
SIM. In our forthcoming paper we will thoroughly analyze the method on more realistic simulations 
and a variety of different planetary systems. 
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Fig. 1. — Solar System Barycenter (SSB) reference frame and the celestial sphere. So is the unit 
vector toward the star with spherical coordinates {a, 5) and ||Ro|| is the distance to the star. 

Fig. 2. — Tangent space at Sq where e^, and are the unit vectors of the local frame. 

Fig. 3. — ASf^ for 150 stars with the largest proper motion from the Hipparcos catalogue. The 
solid lines represent AS^ for 1, 10 and 100 parsecs as a function of Vr Vr. Time span of the mission 
AT = Wyr was assumed. 

Fig. 4. — An^ for 150 stars with the largest proper motion from the Hipparcos catalogue. The 
solid lines represent AII^ for 1, 10 and 100 parsecs as a function of V. Time span of the mission 
AT = lOj/r was assumed. 

Fig. 5. — AS'(^) for 150 stars with the largest proper motion from the Hipparcos catalogue. The 
solid lines represent AS^^^ for 1, 10 and 100 parsecs as a function of V. Time span of the mission 
AT = Wyr was assumed. 

Fig. 6. — Second order effects in the simulated delays for the relative measurements between v 
And and HD 10032 for the baseline vector orientations Bi (a) and B2 (6); (c,d) the same effects 
when the radial velocity of both stars is zero; {e,f) the residuals from the least-squares fit of the 
first order model (69) to the simulated data used in {a,b). 

Fig. 7. — The amplitudes of subsequent harmonic terms for the relative delays Di,D2 (left and 
right panel respectively) corresponding to the planet I (a), H (b) and the parallactic motion (c). 

Fig. 8. — Subsequent steps of the Frequency Decomposition for the simulated relative delay mea- 
surements between v And and HD 10032 corresponding to the baseline vector orientation Bi. Left 
panel contains the residuals after removal of all components from the steps above. Right panel 
contains normalized periodograms of these residuals. 



Fig. 9. — Continuation of Fig. 8 
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Table 1. Target star — v Andromedae (HD 9826, HIP 7513) 



Parameter 


V And 


Right ascension, a (J1991.25) 


0l''36'"47f98 


Declination, 5 (J1991.25) 


41°24'23'/00 


Proper motion in a, iia cos (5 (mas/yr) 


-172.57 


Proper motion in 5, jis (mas/yr) 


-381.01 


Parallax, tt (mas) 


74.25 


Distance, dpc (pc) 


13.47 


Transverse velocity, Vt (km/s) 


26.7 


Radial velocity, Vr (km/s) 


-27.7 


Mass, (Mq) 


1.3 


Orbital elements 


Planet I Planet II 


Semi-major axis, a (AU) 


0.83 2.5 


Semi-major axis, a = -Kam/Mi, (mas) 


0.133 0.813 


Orbital period, P (d) 


241.2 1266.6 


Eccentricity, e 


0.18 0.41 


Epoch of periastron, Tp (JD) 


2450154.9 2451308.7 


Longitude of periastron, a; 


243.°6 247.°7 


Longitude of ascending node, f2 


30.°0 60.°0 


Inclination, i 


45.°0 45.°0 


Mass, m {Mjup) 


2.95 5.98 



Table 2. Reference star — HD 10032 (HIP 7672) 



Parameter 


HD 10032 


Right ascension, a (J1991.25) 


0l''38"48f07 


Declination, 5 (J1991.25) 


40°45'38^'80 


Proper motion in a, ^aCOS(5 (mas/yr) 


-14.70 


Proper motion in 5, fis (mas/yr) 


-2.66 


Parallax, tt (mas) 


8.05 


Distance, dpc (pc) 


124.22 


Transverse velocity, Vr (km/s) 


8.80 


Radial velocity, Vr (km/s) 


-34.00 
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Table 3. SIM orbital elements in SSB reference frame 



Parameter SIM 



Semi-major axis, a (AU) 0.995 

Orbital period, P (d) 362.5 

Eccentricity, e 0.015 

Epoch of periastron, Tp (JD) . . 2451519.44 

Longitude of periastron, 6<j 74.°67 

Longitude of ascending node, O, 0.°005 

Inclination, i 23.° 45 



Table 4. Dominant planetary terms from Frequency Decomposition 



Parameter Planet I Planet II 



/ (1/d) .. 1/241.35 1/1265.87 

C'l (m) ... -0.743 X IQ-^^ 0.282 x 10"^ 

Si (m) ... 0.586 X 10"^ 0.790 x lO"" 

(m) ... -0.481 X 10-^ 0.756 x 10"^ 

(m) ... 0.146 X 10-8 0.329 x 10"^ 

df (m) ... 0.198 X 10-9 -0.301 x 10"^ 

Sf (m) ... -0.683 X 10"^ 0.466 x 10"^ 

C| (m) ... 0.492 X 10"^ -0.624 x 10"^ 

5| (m) ... -0.148 X lO-'' -0.202 x 10"^ 



-23- 



Table 5. Orbital elements derived from FD parameters 



Parameter 


Planet I 


Planet II 


Semi-major axis, a (mas) 


0.130 


0.733 


Orbital period, P (d) 


241.35 


1265.87 


Eccentricity, e 


0.22 


0.39 


Epoch of periastron, Tp (JD) . . 


2450164.79 


2451301.97 


Longitude of periastron, lo 


255.°05 


243.° 15 


Longitude of ascending node, CI 


31.°09 


62.°95 


Inclination, i 


44.° 50 


43.° 13 



Table 6. Orbital elements from standard model 



Parameter 


Planet I 


Planet II 


Semi-major axis, a (AU) 


0.132 


0.813 


Orbital period, P (d) 


241.21 


1265.65 


Eccentricity, e 


0.17 


0.41 


Epoch of periastron, Tp (JD) . . 


2450152.63 


2451306.95 


Longitude of periastron, lo 


240.°34 


247.° 88 


Longitude of ascending node, CI 


29.°65 


59.°96 


Inclination, i 


44.°91 


45.°04 
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Fig. 2.— 
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